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Finite-temperature (T) properties of a Kitaev model defined on a honeycomb lattice are investigated by a 
quantum Monte Carlo simulation, from the viewpoint of fractionalization of quantum 5 = 1/2 spins into two 
types of Majorana fermions, itinerant and localized. In this system, the entropy is released successively at 
two well-separated T scales, as a clear indication of the thermal fractionalization. We show that the high-T 
crossover, which is driven by itinerant Majorana fermions, is closely related with the development of nearest- 
neighbor spin correlations. On the other hand, the low-T crossover originates in thermal fluctuations of fluxes 
composed of localized Majorana fermions, by which the spectrum of itinerant Majorana fermions is significantly 
disturbed. As a consequence, in the intermediate-T range between the two crossovers, the system exhibits T- 
linear behavior in the specific heat and coherent transport of Majorana fermions, which are unexpected for 
the Dirac semimetallic spectrum in the low-T 1 limit. We also show that the flux fluctuations tend to open an 
energy gap in the Majorana spectrum near the gapless-gapped phase boundary. Our results indicate that the 
fractionalization is experimentally observable in the specific heat, spin correlations, and transport properties. 

PACS numbers: 75.10.Kt, 75.70.Tj, 75.10.Jm 


I. INTRODUCTION 

The fractionalization of electrons in solids is one of the cen¬ 
tral topics in modern condensed matter physics. A prototyp¬ 
ical example is found in one-dimensional strongly correlated 
electron systems: charge and spin degrees of freedom in an 
electron behave as independent particles, which are termed 
holon and spinon, respectively [JJ]. A different form of frac¬ 
tionalization is also anticipated in insulating magnets with ge¬ 
ometrical frustration. For instance, the existence of the ele¬ 
mentary excitation carrying a half of spin, named spinon, is 
predicted in a quantum spin liquid (QSL) @-01. and emer¬ 
gence of magnetic monopoles is suggested in spin ice sys¬ 
tems H- Another fractionalization was pointed out in heavy 
fermion systems as well. The half residual entropy in the two- 
channel impurity Kondo system is understood from the frac¬ 
tionalization of S' = 1/2 impurity spin into two Majorana 
fermions @], 

A quantum spin model, called the Kitaev model, has re¬ 
cently attracted considerable attention in broad areas of re¬ 
search, not only condensed matter physics but also statistical 
physics and quantum information OQ - This model is com¬ 
posed of S = 1/2 spins with bond-dependent interactions on a 
honeycomb lattice. Such peculiar interactions were suggested 
to be realized in the systems with strong spin-orbit coupling, 
such as iridium oxides |H]. The most striking feature of this 
model is that it is exactly solvable due to the existence of Z 2 
conserved quantity on each hexagon, termed flux. The ground 
state dictates both gapless and gapped QSL phases depend¬ 
ing on the exchange coupling constants. The exact solution 
is provided by representing 5 = 1/2 spins by two types of 
Majorana fermions: one is localized and composes the fluxes, 
and the other forms itinerant bands mn. The latter itiner¬ 
ant Majorana fermions determine the excitation spectrum in 
the QSLs. Thus, the fractionalization of spins into Majorana 
fermions is not just a mathematical tool but physically impor¬ 


tant in the Kitaev model. 

A natural question arising here is how high-temperature (T) 
paramagnetic spins are fractionalized into Majorana fermions 
when cooling the system. Considering the fact that the spin- 
charge separation in the one-dimensional electron systems 
plays a key role in comparison with experiments, it is cru¬ 
cial to elucidate the thermal fractionalization for experimental 
exploration of the QSL physics. The thermodynamic proper¬ 
ties in the Kitaev model and its extensions have been studied, 
mainly for explaining the magnetism in iridium oxides 1IT2I — 
El. but the signature of fractionalization at finite T was not 
addressed in most of the previous studies. Among them, how¬ 
ever, the authors pointed out the significance of fractional¬ 
ization in a peculiar phase transition at finite T in a three- 
dimensional extension of the Kitaev model l T5lfl6ll : the phase 
transition is governed by thermal excitations of localized Ma¬ 
jorana fermions. Nevertheless, the relevance of fractional¬ 
ization remains unclear, in particular, to the experimentally- 
observable quantities. 

In this paper, we investigate the effect of fractionalization of 
quantum spins on the finite-T properties of the Kitaev model 
on a honeycomb lattice by applying the unbiased quantum 
Monte Carlo (QMC) method. In this model, the two Majo¬ 
rana fermions, itinerant and localized, release their entropy 
successively at two well-separated T scales. We elucidate that 
each crossover has an impact on experimental observables: 
the high-T one, driven by itinerant Majorana fermions, corre¬ 
sponds to the development of spin correlations between neigh¬ 
boring sites, while the low-T one, originating from thermal 
fluctuations of localized Majorana fermions, is accompanied 
by a sizable change in the excitation spectrum of itinerant Ma¬ 
jorana fermions. This leads to apparent T-linear behavior of 
the specific heat and coherent transport of Majorana fermions 
in the intermediate-T state between the two crossovers, in 
contrast to the Dirac semimetallic behavior and T 2 specific 
heat anticipated in the low-T limit. Moreover, we show that 
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the thermal excitation of fluxes tends to open a gap at finite T 
near the gapless-gapped phase boundary. 

The paper is structured as follows. In Sec. El we introduce 
the Kitaev model on a honeycomb lattice, and briefly review 
the ground-state properties. In Sec. [Till we present the numer¬ 
ical method to analyze the finite-T properties in the Kitaev 
model. The definitions of physical quantities are also given 
in this section. The numerical results are shown in Sec. m 
We present the T dependences of the specific heat and the en¬ 
tropy, and summarize the two crossovers in the phase diagram 
in Sec. II V Al We discuss the origins of the crossovers by cal¬ 
culating the spin correlation in Sec. II V Bl and the flux density 
in Sec. II V Cl We also compute the density of states (DOS) of 
the itinerant Majorana fermions in Sec. HV Dl We discuss the 
peculiar T dependence of the specific heat in the intermediate- 
T region in Sec. IIV El In Sec. IIV FI we evaluate the optical 
conductivity and the Drude weight of the itinerant Majorana 
fermions at finite T. The effect of thermal fluctuations near 
the gapless-gapped boundary is discussed in Sec. IIV Gl Fi¬ 
nally, Sec.[V]is devoted to the summary. 

H. MODEL 

The Kitaev model is composed of S' = 1/2 spins defined 
on a honeycomb lattice, whose Hamiltonian is given by 0 

H = -J x ]T &jo% ~J y J2 ° V A J* E (D 

{jk)x (ffe)v OUz 

where erj is the /(= x,y,z ) component of the Pauli matrix 
representing an S = 1/2 spin at site j. Corresponding to 
three inequivalent bonds on the honeycomb lattice, named x, 
y, and z bonds, the sum over (jk)i is taken over the nearest 
neighbor (NN) sites on the l bonds. 

The ground state of the model in Eq. CD was exactly solved 
by introducing Majorana fermions 0. The ground state has 
gapped and gapless excitations depending on the exchange 
constants, J x , J y , and J z 01 (see the inset of Fig. [6]t. The spin 
correlations are extremely short-ranged, i.e., nonzero only for 
the NN pairs for all the parameters, which indicates that both 
gapped and gapless ground states are QSLs lfl8i - l20ll . The 
model does not exhibit any phase transition at finite T al¬ 
though a three-dimensional variant does 10- 

Hereafter, we describe the anisotropy of the exchange con¬ 
stants by the parameter a as J x = J y = a/3 and J z = 
1 — 2a/3 (J x + J y + Jz = 1) as shown in the inset of 
Fig. [6] Along this cut in the ground-state phase diagram, the 
gapped-gapless phase boundary in the ground state is located 
at a = 3/4. 


HI. METHOD 

An exact solution for the ground state of the Kitaev model 
is formulated by the Jordan-Wigner transformation along the 
chains consisting of the x and y bonds 00- The fermions 
introduced by the transformation can be represented by two 


Majorana fermions Cj and c :l at each site j. Using these Ma¬ 
jorana fermions, the Kitaev model is rewritten as 

kL — iJx ^ ' CjCk 1 J y / ' CjCfc iJz ^ ~ CjrCjCki ( 2 ) 

{jk) x {jk)y (jk) z 

where the sum over (jfc) is taken for the NN sites with j < k. 
The operator r/ r = iCjCk is defined on each z bond (r is the 
bond index). This is regarded as a classical variable taking ±1 
because of [H , rj r ] = 0 and rj x = 1 for all r. 

By using this Majorana representation, we carry out the 
QMC simulation at finite T. For the model in Eq. ([J), the 
partition function Z is written in the form 

Z = Tr M Tr {ci} e- 0n = E (3) 

{r) r }=±l 

where Tr{ r/r . j and Tr/ Cj \ are the traces for localized and itin¬ 
erant Majorana fermions, respectively; /3 = 1/T is the in¬ 
verse temperature (we set the Boltzmann constant k]>, = 1). 
Here, Ff({r] r }) is the free energy of the Majorana fermion 
system for a fixed configuration {ry r }, which is easily calcu¬ 
lated by the exact diagonalization. We perform the Markov 
chain Monte Carlo (MC) simulation for sampling the con¬ 
figurations of {rj r } so as to reproduce the thermal distribu¬ 
tion of e~P F f J4 '7’*F) [16]. In the present calculations, we 
performed the QMC simulation hybridized with the paral¬ 
lel tempering technique with 16 replicas tf2~ij| . We spent the 
10,000 MC steps for thermalization and 40,000 MC steps for 
measurement in up to an L = 12 cluster, which contains 
N = 2 x L 2 = 288 sites. 

In the QMC simulation, we calculate the specific heat as 



where Ef({r/ r }) is the energy of the itinerant Majorana 
fermion system for a given configuration {rj r }. From the T 
dependence of the specific heat, we obtain the entropy per site 
as 

S = \n2~J m dT'C v /T', (5) 

where T m is chosen to be T m = 10 (3> J x + J y + J z = 1). 

We also calculate the equal-time spin correlations. In the 
Kitaev model, as (cjcr^,) ^ 0 forNN bonds lfl8ll . we measure 
the NN spin correlations by 

< 6 > 

(jk)i 

which are given by each term in Eq. (J2]» in terms of the Ma¬ 
jorana fermions. In addition, we compute the thermal average 
of the flux density as 

W =|fETO’ C7) 

P 
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where W p is composed of rj T included in the hexagonal pla- 
quette p: W p = UrepVr- 

In addition to the above thermodynamic quantities, we 
calculate the dynamical quantities for itinerant Majorana 
fermions. The DOS of the itinerant Majorana fermions with a 
given configuration of rj r is defined by 

D(u},{r] r }) = Y S (“ - E n {{r] r })), ( 8 ) 

n 

where e n is the one-particle energy of the fermion /,, which 
is introduced so as to diagonalize the Hamiltonian as 

KiM) = Y £ n{M) (flfn ~ ^ ■ ( 9 ) 

n ' ' 

We calculate the thermal averages of the DOS, D(w) = 
{D(oj, {?7r})} by the QMC simulation. Note that D(oj) do not 
contain the T dependence of the Fermi distribution function: 
we take into account the effect of thermal fluctuations only on 

Vr- 

Moreover, we compute the optical conductivity of itinerant 
Majorana fermions. For this purpose, we introduce the Fourier 
transform of the Hamiltonian as 

H = Y, c k^kCk = Y £ -k (flkf nk ~ 2 ) ’ 

k n:£nk> 0 k ' ' 

( 10 ) 

where Ck is a set of the Fourier transforms of Cj and the Lx L 
cluster is regarded as a unit cell. The Bloch Hamiltonian Hk is 
diagonalized by introducing a set of fermions f n k belonging 
to the n-th band with the energy Then, the conductivity 
tensor is calculated by 

1 r°° W3 

a^iu) = T dte^^ d\(M-iX)Mt)), ( 11 ) 
L Jo Jo 

where 6 is an infinitesimal positive number, 0(t) = 

anc | curren t operator is defined as J= 

Yjknn’ fknf k n'{ukn\dHk/dk p ,\u kn ') with the eigenstate 
| Ukn) of Hk- To extract the contribution to coherent transport, 
we also obtain the Drude weight via the sum rule. Specifically, 
we compute the Drude weight along the x direction by 

i i r°° 

D * = wY ~ - / (!2) 

<«)i n Jo 

where the summation J2uj)' I s taken only for the NN x bonds 
on the boundary of a finite-size cluster which we regard as a 
unit cell in the calculations. 

IV. RESULTS 
A. Specific heat and entropy 

Figures |Xt a )-tIl d) show the QMC data for the specific heat 
C v [Eq. 01 as a function of T for several values of the 


anisotropy parameter a. For all cases, the specific heat ex¬ 
hibits two peaks; both are almost system-size independent, in¬ 
dicating two crossovers. We hereafter term the low- and high- 
T crossover temperatures as Tl and Th, respectively. Fig¬ 
ures [D e)0h) show the entropy per site obtained by Eq. 0. 
The entropy rapidly decreases with decreasing T in the vicin¬ 
ity of Tl and Th corresponding to the two peaks of the specific 
heat. A half of the entropy is released successively in each 
crossover; consequently, the entropy becomes ~ | In 2 per 
site in the region between Tl and Th- The plateau-like behav¬ 
ior of the entropy in this region becomes clearer for smaller a, 
i.e., larger anisotropy of the exchange constants. 

Figure [2] shows the contour map of the entropy on the a-T 
plane. The high-T crossover temperature Th is almost inde¬ 
pendent of a. The origin will be discussed in the next section 
II V Bl On the other hand, the low-T crossover temperature Tl 
strongly depends on a\ this will be discussed in Sec. IIV Cl 
There are three regions separated by the two crossovers Tl 
and Th- In the following sections, we clarify the differences 
between these regions and their influences on the observable 
quantities. 


B. High-T crossover: spin correlation 

Let us first discuss what takes place in the high-T crossover 
at Th- The itinerant Majorana fermions c :/ form a band whose 
width is Wb = 2 (J x + J y + J z ) — A = 2 — A, where A is the 
excitation gap in the gapped phase. Suppose that the system 
is in the gapless region and the DOS is constant ~ 1 /Wb, the 
specific heat originating from the itinerant Majorana fermions 
takes maximum at T ~ 0.511. This value well coincides with 
Th in a wide region of a, even in the gapped region for a < 
0.75, as shown by the dashed-dotted line in Fig. [2] The result 
clearly indicates that the high-T crossover originates in the 
itinerant Majorana fermions. 

We find that the crossover at Th is closely related with the 
development of NN spin correlations given inEq. 0, which is 
observable in experiments. The T dependences of S u are pre¬ 
sented in Figs.0i)01) for the same set of a as in Figs. 0a)- 
[U h). In the high-T limit, S l j k is given by the high-T expansion 
as Tr[crjcr[e _/3w ]/Tre _ ^ w ~ —/3Tr[<r \o l k W[ = f3Ji. Our 
QMC data obey this Curie behavior, indicated by the dashed- 
dotted curves in Figs. IHi)-[Ill). In the crossover region near 
Th, however, the spin correlations show deviations from the 
Curie behavior, and quickly saturate to the values that are ana¬ 
lytically obtained for the ground state (horizontal dashed lines 
in the figures) |18|. Hence, the high-T crossover by the itin¬ 
erant Majorana fermions corresponds to physically important 
behavior in this quantum spin system: the growth of the NN 
spin correlations. We note that the spin correlations also show 
slight changes in the low-T crossover at Tl- This behavior is 
discussed in Sec. IIV Gl 
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FIG. 1: (color online), (a)-(d) T dependences of the specific heat at (a) a = 1.0, (b) a = 0.8, (c) a = 0.75, and (d) a = 0.5 in the 
several clusters with 2 x L 2 spins. Here, we define the anisotropy parameter a by taking J x = J y = a/ 3 and J z = 1 — 2a/3. (e)-(h) T 
dependences of the entropy per site, S , and the thermal average of the density of the flux W p , W. (i)-(l) T dependences of the equal-time spin 
correlations, S 11 ', S p = (S xx + S yy )/ 2. The horizontal dashed lines represent the values at T = 0 which are calculated analytically d, and 
the dashed-dotted curves represent the high-T Curie behaviors S 11 ~ Ji /T. 
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FIG. 2: (color online). Contour map of the entropy per site, S/ In 2, 
on a plane of T and a. The dashed line represents crossover temper¬ 
ature obtained by the perturbation theory in the limit of J z J x , J y 
(a <C 1). The dashed-dotted line represents the crossover tempera¬ 
ture obtained by assuming the constant DOS. NNC stands for a NN 
correlation. 


C. Low-'/' crossover: flux density 


Next, we discuss what occurs in the low-T crossover. The 
entropy release near Tl originates from the localized Majo- 
rana fermions Cj or q, r . This is confirmed by calculating the 
T dependence of the flux density, W, in Eq. 0, as shown 
in Figs.|TJe)0h). The results show that W rapidly decreases 
from 1 with increasing T in the vicinity of Tl. Hence, the 
crossover at Tl is due to the thermal fluctuation of fluxes. 

This is further confirmed by considering the toric code limit 
corresponding to J x ,J y <C J z (a <C 1). In this limit, 
the Kitaev model is reduced to the effective model T-L e s = 


Jeff J2 P Wp> where J eS = J 2 Jy/{16 J 3 ) 0]. Since this ef¬ 
fective model describes free Ising spins in the magnetic field 
Jeff, the specific heat is of Schottky-type, which takes a max¬ 
imum at Tl/J e ff ~ 0.833. The asymptotic behavior is shown 
by the dashed line in Fig. [2] The agreement between this line 
and Tl further supports that the low-T crossover is caused by 
the thermal disturbance of fluxes. 

We also note that the agreement of the asymptotic behav¬ 
ior is consistent with the absence of phase transition in this 
two-dimensional system. This is in contrast to the three- 
dimensional case; there are local constraints for W p in the Ki¬ 
taev model defined on a three-dimensional hYperhoneycomb 
lattice, leading to a finite-T phase transition ltl6i l. On the other 
hand, there is no constraint for W p in the two-dimensional 
case, which results in the absence of the phase transition for 
T > 0. 

To summarize the above results, the three regions in the 
phase diagram depicted in Fig.[2]are characterized as follows. 
The high-T region for T > Th is a conventional paramagnetic 
state, where the NN spin correlations obey the Curie behavior. 
On the other hand, in the low-T region for T < Tl, the NN 
spin correlations are saturated to the T = 0 values, and fur¬ 
thermore, the fluxes are also aligned. Thus, the system below 
Tl behaves similar to the ground state QSL. In the region for 
Tl < T < Th, a peculiar intermediate state appears: the NN 
spin correlations are well developed, whereas the fluxes are 
thermally disordered. In the following sections, we discuss 
the nature of this intermediate state. 


D. Density of states of itinerant Majorana fermions 

Since the Z 2 variables r/ r couple with the itinerant Ma¬ 
jorana fermions, we expect that the enhanced fluctuations 
of fluxes above Tl affect the nature of itinerant Majorana 
fermions considerably. In order to elucidate such behavior, we 
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FIG. 3: (color online). The DOS of Majorana fermions at (a) a = 1.0, (b) a = 0.9, (c) a = 0.8, and (d) a = 0.75. Except the results at 
T = 0 and T = oo, the DOS are calculated by QMC for the 10 x 10 superlattice of the L = 12 cluster. 


gion, the low-7' specific heat is expected to be proportional to 
T 2 because of the Dirac semimetallic dispersion for aligned 
fluxes. However, C v calculated by assuming all rj r = +1 
largely deviates from our QMC data in the calculated T range, 
as shown for a = 1.2 in Fig. [4] This indicates that the asymp¬ 
totic T 2 behavior will be limited only in the extremely low-T 
region, much lower than Tl. 

Instead, in a wide range of Tl < T < Th, we find that C v 
well scales to oc T, which originates from the “metallic” DOS 
caused by thermally fluctuating fluxes above Tl. Indeed, the 
overall behavior including T > Th is well explained by the 
result for completely random {r] r }, as shown in Fig. [4] Thus, 
as a consequence of the thermal fractionalization of quantum 
spins, we find the apparent T-linear behavior, not T 2 , in the 
intermediate-T region where the NN spin correlations are well 
developed. 


F. Optical conductivity and Drude weight 

The significant change in the DOS of the itinerant Majo¬ 
rana fermions will affect transport properties as well. We here 
show it by computing the optical conductivity of itinerant Ma¬ 
jorana fermions [Eq. (ITTb ]. Here, we consider the longitudinal 
component along the x direction (//. = v x). The calcula¬ 
tions were done for the lxl supercell of the L = 12 cluster. 

Figure 0a) shows the results of <j xx (u>) at several T for 
a = 1.0. The incoherent component at finite ui increases 
with decreasing T below Th- To extract the contribution to 
coherent transport, we calculate the Drude weight D x of the 
itinerant Majorana fermions by using the sum rule in Eq. (IT2l) . 
Figure 0b) shows the T dependence of D x . While decreas¬ 
ing T, the Drude weight gradually increases below Th, and 
sharply decreases to zero below Tl after showing a peak near 
Tl. The result suggests that the transport quantities, such as 
the thermal conductivity, have sizable values between the two 
crossovers. 


G. Gapless-gapped phase boundary 

Finally, we discuss the effect of thermal fluctuations near 
the phase boundary between the gapless and gapped phases. 
Whereas the coherent transport and the T-linear behavior are 
observed widely in the region where the ground state is gap- 
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FIG. 4: (color online). T dependence of the specific heat C v at a = 
1.2 in the L = 12 cluster. For comparison, the results calculated by 
fixing all T] r to +1 and by assuming random {r/ r } are shown by the 
solid and dashed curves, respectively. The log plot of C v /T is also 
shown in the inset. 


calculate the DOS of itinerant Majorana fermions. The calcu¬ 
lations were done for the 10 x 10 supercell, where the L = 12 
cluster obtained by the MC simulation is regarded as a unit 
cell. The calculations at T = 0 (T = oo) are performed for 
a L = 6, 000 (X = 60) cluster. In the calculation at T = oo, 
we take a simple average over 10,000 random configurations 
of {rjr}. 

Figure 0a) shows the result for the isotropic case a = 1.0 
(J x = J y = ./,). The QMC data are shown near Tl, to¬ 
gether with the results at T = 0 and T = oo. In this gap¬ 
less case, at T = 0, the DOS shows semimetallic behav¬ 
ior D(u}) = (D(co, {rj r })) oc oj for small u>, reflecting the 
Dirac dispersion. While increasing T above Tl, however, the 
semimetallic dip of DOS is filled rapidly, leading to “metal¬ 
lic” behavior, D[u) = 0) ^ 0. The result clearly indicates that 
the thermal fluctuations of fluxes near Tl significantly affect 
the low-energy spectrum of itinerant Majorana fermions. 


E. Peculiar T dependence of the specific heat 

The significant change in the DOS for T > Tl affects the 
T dependence of the specific heat C v . In the gapless QSL re¬ 
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FIG. 5: (color online), (a) The optical conductivity at a = 1.0 on the 
L = 12 cluster at several T. (b) T dependence of the Drude weight 
of itinerant Majorana fermions at a = 1.0. 


less, they are disturbed in the vicinity of the phase boundary at 
a = 0.75. In this region, thermal fluctuations in the Z 2 vari¬ 
ables { Tj r } bring about different behavior in the low-energy 
part of D(u>) = (D(uj, {f] r }))- Figures |2b), 0c), and0d) 
show the DOS of the itinerant Majorana fermions at a = 0.9, 
a = 0.8, and a = 0.75, respectively. At a = 0.8 and 
a = 0.75, the system develops an energy gap with increas¬ 
ing T in the vicinity of Tl, in sharp contrast to the gap filling 
at a = 1.0 and a = 0.9. The results indicate that there is an 
intermediate region where the thermal fluctuation of rj r gaps 
out the low-energy excitation of itinerant Majorana fermions. 

The intermediate region is identified by calculating the 
magnitudes of the gaps at T = 0 and T = oo, as presented in 
Fig. [6] The schematic phase diagram determined by the DOS 
at T = oo is presented in the inset. Remarkably, the gapped- 
gapless boundary is similar to that in the dynamical phase di¬ 
agram |j22i|, suggesting a relation between thermal and quan¬ 
tum fluctuations. We also note that the boundary is similar to 
the result for the full flux state ll23tl . 

The modification of the boundary at finite T implies that ef¬ 
fective exchange couplings are renormalized in an anisotropic 
way by the thermal fluctuation of the Z 2 variables {//,.}. In¬ 
deed, the anisotropy of spin correlations is slightly enhanced 
near Tl while increasing T, as shown in Figs. |T]j) and|T]k). 
Thus, the slight change in the spin anisotropy is interpreted 
as a consequence of the change of the excitation gap in the 
itinerant Majorana fermions fractionalized from the quantum 
spins. 


V. SUMMARY 

In summary, we have investigated the thermal fractionaliza- 
tion of quantum spins into Majorana fermions in the Kitaev 
model by using the QMC simulation. We clarified that the 
fractionalization appears as two crossovers, both of which are 
physically observable in the thermodynamics. The higher-T 
crossover is identified by the development of NN spin cor¬ 
relations, which will be observed in, e.g., neutron scattering 
experiments. In between the crossovers, the Drude weight 
of itinerant Majorana fermions takes a sizable value, which 
might be observed by the thermal conductivity. Meanwhile, 
the low-T crossover leads to a peculiar T-linear behavior in 
the specific heat above the crossover temperature. We also 
showed that the thermal fractionalization affects the gapped- 
gapless phase boundary by renormalizing the spin anisotropy. 
The present results complete how the fractionalization of 
quantum spins into Majorana fermions occurs while changing 
temperature in the ideal Kitaev model. This provides a use¬ 
ful reference to the experimental exploration of QSLs in, e.g ., 
iridium oxides 124142811 and ruthenium compounds I29l432ll . 
where the dominant interaction is expected to be of Kitaev 
type. 
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FIG. 6: (color online). The excitation gap for the Majorana fermions 
at T = 0 (blue solid line) and T = 00 (red symbols) as a function 
of a. The inset indicates the gapped-gapless boundaries on the plane 
of J x + J y + Jz = 1. The blue solid lines represent the phase 
boundaries in the ground state, while the red dashed lines represent 
the boundaries obtained from the DOS at T = 00 . See the text for 
details. 
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